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ABSTRACT 

The objective of this paper is to find the best-fit probability model for annual maximum rainfall data for Botswana 
for the period January 1901 to December 2012. The Kwiatkowski-Philips-Schmidt-Shin test of stationarity on the data 
shows that the time series is stationary. Based on visual inspection of the generalized quantile plot, three extreme value 
probability distributions are considered, belonging to the Gumbel and Frechet maximum domains of attraction: the 
gamma, lognormal and BurrXII. The maximum likelihood method is used to estimate the parameters of each distribution 
under study. The empirical CDF plots and Q-Q plots of the data show that the three models are highly competitive in terms 
of goodness-of-fit. Formal model assessment criteria, namely, the Anderson-Darling test and the Bayesian Information 
Criterion agree on the ranking of the three models: the lognormal distribution gives the best fit to the data followed closely 
by the gamma whilst the BurrXII distribution comes a distant third. 
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INTRODUCTION 

Precipitation and temperature extremes are considered to be the most important climate events and have been 
extensively explored over the past several decades, according to [1] and [2]. Water is one of the most important natural 
resources with diverse uses in the agricultural, domestic and industrial sectors. Low or excessive rainfall can lead to water 
scarcity or floods, which may cause great damage to human health, the environment, and property, leading to massive 
social and economic losses. While nearly a third of the people in Africa already live in the drought-prone areas, it is 
estimated that climate change will add up to over 80 million people at risk of hunger by 2080 [3-4], The debilitating effects 
of climate change have not spared Botswana, which is characterized by a semi-arid climate with an erratic and low rainfall 
(450 mm per annum), high evaporation (more than 2000 mm per annum) and limited fresh water supplies. According to 
[5], typically, between 1965 and 2006 the country has suffered from several droughts and flood-related disasters which 
have affected 1.5 million people and caused economic losses of over US$ 3 million to the country. Most of the water used 
for domestic, industrial and agricultural purposes in Botswana is from dams. Consequently, water demand far exceeds 
water supply as the dams are drying up due to erratic and low rainfall coupled with high evaporation rates across the 
country. This situation is a major threat to the livelihoods of the country’s population, especially the rural stratum. 
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According to [6], about 70% of the rural households derive their livelihoods from agriculture through subsistence farming 
and crop production is mainly based on rain-fed farming. 

The main motivation for statistical modelling of weather extremes is reliability: the ability of a system or 
component to withstand the weather conditions for the whole of its working life. For instance, for civil engineers to 
determine how high the wall of a reservoir dam should be, they need to estimate what the heaviest (maximum) rainfall 
amount will be over some specified future time period. This can only be achieved by using historical data and fitting an 
appropriate statistical model for them. That is, a model for risk is developed by selecting a probability distribution that 
gives the best fit to the data at hand. Statistical models for extreme values are based on an asymptotic theory called extreme 
value theory. Extreme value theory provides the statistical framework to make inferences about the probability of very rare 
and extreme events [7]. The theory is based on the probability distribution of the unusually low or large values of a time 
series. 

The application of extreme value theory to the analysis and modelling of extreme rainfall from different regions of 
the world bears testimony to the importance of determining the frequency of extreme rainfall events in different parts of the 
world. For instance, [8] determine the best-fit models for the frequency analysis of extreme monthly rainfall events in 
Bangladesh and estimate extreme quintiles for several periods of time; [7] provide the first application of extreme value 
distributions to model minimum annual rainfall in Zimbabwe;[4] model a long-term annual rainfall data from 1961 to 2003 
at 11 synoptic stations across Botswana using extreme value distributions.[9] apply extreme value theory to model rainfall 
data from South Korea; [10] use extreme value theory to model rainfall data from Europe and the USA; [11] apply extreme 
value theory to Greece’s rainfall data; while[12] model rainfall data from Italy using extreme value distributions. 

Though modelling of extreme rainfall has been shown to be important, not much research in this area has been 
done in Botswana. Perhaps the most significant research so far known to the author is work by [4], which involves a 
regional frequency analysis of annual rainfall data of covering 40 years at 11 synoptic stations across Botswana using 
extreme value distributions. The results of the study show that all the stations behaved homogeneously and followed the 
generalized extreme value (GEV) distribution. However, as pointed out by [13], rainfall should be monitored overtime as 
significant changes in the rainfall may lead to floods or drought. Knowledge of the distribution, intensity, and duration of 
precipitation extremes is important to stakeholders like climate change researchers, businesses, governments and members 
of the general public as forecasting of future events may help them adequately prepare for any extreme events. 

The main objective of the present study is to find the best-fit model for determining the frequency of extreme 
rainfall events in Botswana by using historical annual maximum rainfall data from 1901 to 2012. The other objective seeks 
to determine whether there has been any monotonic change (trend) in the annual maximum average rainfall in Botswana in 
the period 1901-2012. We compare the performance of three different distributions: the gamma and lognormal selected 
from the Gumbel and, the BurrXII from the Frechet maximum domains of attraction of normalized maxima.. The 
remainder of the paper is organized as follows. The second section provides a description of the annual maximum rainfall 
data in Botswana from 1901 to 2012 and the extreme value methodology employed in the study. The third section presents 
the results and a discussion of the study. The fourth section contains the conclusions of the study. 
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DATA AND RESEARCH METHODOLOGY 

Data Description 

The data used for this study consist of annual maximum rainfall for Botswana extracted from the mean historical 
monthly rainfall for Botswana during the time period 1901-2012. The dataset was produced by the Climatic Research Unit 
(CRU) of the University of East Anglia (UEA), United Kingdom, and is available at the World Bank website 
http://sdwebx.worldbank.org/climateportal/index.cfm?page=downscaled data download&menu=historical . The choice of 
the time period is based on the quality and availability of data. The data produced by the CRU are of very high quality with 
no missing values. 

Theoretical Foundation of Extreme Value Theory 

The basis for extreme value theory is the Fisher-Tippet-Gnedenko Theorem (Extremal Types Theorem), first 
discovered by [14] and later proved in full generality by [15]. The theorem postulates the asymptotic distribution of 
extreme order statistics. It states that if , X 2 ..., X n is a sample of independent and identically distributed (iid) random 

variables with a common distribution function F, and X {1) , X^y.., X (n) is the increasingly ordered sample, and 

X (n) =max(Xp X 2 ..., X n ) is normalized with constants a n > 0 (for all n) and b n such that 


lining P 


( X{n) < z) = lim n -»oo(E(a n z + b n )} n = G f (z) for all z £ 

V a n / s 


( 2 . 1 ) 


where Gr(z) is a non-degenerate distribution function, then GXz) is called an extreme value distribution function. It is 


given by Equation (2.2). 
G^(z) = 


jexp(—(1 + fz) f), 1 

(. exp(— exp(z)), 


+ fz > 0 if cf & 0 
z £ M if = 0 


( 2 . 2 ) 


exp(— exp(z)), 

where f is called the extreme value index (EVI), a parameter that determines the heaviness of the tail of the distribution. 
And, the only possible types for Gc(z) (limiting distributions for the linearly normalized variable Z — —^——) as defined 
in Equation (2.1) are the Gumbel (which has an infinite right endpoint with the tail that decays much faster than the tail of 
Frechet) with distribution function G^(z) = exp exp (~)) < z E R ; Frechet (with an infinite right endpoint) having the 


distribution function Gc (z) = 


0 , 


z < b 


exp 


(-(=?)■'). ■>» 


, and the Weibull (short-tailed with a finite right endpoint) with 


G^(z) - 


exp 


(-[-ran- ■<». 

z > b 


1 , 


Important to note is that all the three types of distributions have a location and scale parameters, b and a, with a > 0, 
respectively. In addition, the Frechet and Weibull have a positive shape parameter OC , which gives the EVI through the 
relations f = Voc an ^ f = _ Voc f° r the former and latter, respectively. 
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A reformulation of these three limiting distributions for the normalized maximum M n involves combining them 

to form a single family called the generalized extreme value (GEV) distribution with a cumulative distribution function 
given by equation (2.3). 



{ z; i + f (£_») > 0}, ? *0 
z £ M, f = 0 


(2.3) 


with — oo < /r < oo and a > 0 and — oo < f < oo . The p, o and f are the locations, scale and shape parameters, 
respectively, with f determining the rate of tail decay. The GEV distribution incorporates light-tailed Gumbel for f -» 0, 
heavy-tailed Frechet forf > 0, and short-tailed Weibull for f < 0. For details on introduction to statistical modelling of 
extremes, the reader is referred to [16]. 


Selecting the Maximum Domain of Attraction 

A distribution function F is said to be in the domain of attraction for maxima of G%, written F £ D(G ^), whenever 
condition (2.1) holds. And, as pointed out by [17], testing whether F £ D(G^), for a certain £, is a crucial topic when an 
application of extreme value procedures is to be considered. However, in practice, the analyst has no information on the 
sign of the extreme value index, f, a priori. A brief overview of some testing procedures for the extreme value condition 
and for the statistical choice of the tail are discussed in [17]. However, such procedures are not easy to implement and, in 
practical situations, exploratory graphical tools like the quantile (Q-Q) plot and the mean excess plot are used to determine 
the tail behavior of the underlying distribution of sample data. 

In this paper, for ease of implementation, we use a graphical tool called the generalized quantile plot defined by 
[18]. Given a sample of iid random variables X lt X,..., X n , with X l{) , X (2) ..., X jn) the corresponding increasingly 

ordered sample, the generalized quantile plot is based on the generalized Hill estimator for E, , which is a generalisation of 
the ordinary estimator of the slope of the k largest points of the Pareto Q-Q plot. From [20], the generalized Q-Q plot is a 
graph of the points (log (~~) < log ( X n _ kn H kn ) ) , fc = 1,2,..., n — 1, where X n _ kn H kn , called the UH scores are usually 
written as UHj n — X n _j n H jn , j = 1,2, ...,n — 1, and// k n are the Hill estimates given byH kn — 1/k Yg=i logX n _ ]+ln — 
log UH n _ kn - 

Clearly, as pointed out by [18], for small k the generalized Q-Q plot becomes ultimately linear. To be noted also is 
that the generalized Hill estimator discussed in [18] was proposed earlier as an estimator for the extreme value index 
fusing the slope of the generalized Q-Q plot in the case f £ M. For these reasons, the slope of the generalized Q-Q plot can 
be used as a simple practical tool to check for the maximum domain of attraction condition for all the three types of 
extreme value distributions: the slope of the generalized Q-Q plot in the largest observations is positive if f > 0, indicating 
that the underlying distribution is heavy-tailed (in the Frechet domain); is fairly constant if f = 0, indicating that the 
underlying distribution is medium-tailed (in the Gumbel domain) and is negative if f <0, indicating that the underlying 
distribution is short-tailed with a finite right endpoint (in the Weibull domain). 
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Testing for Stationarity 

The statistical basis for parameter estimation and time series forecasting depends on the series being second-order 
(covariance) stationary: data not having any upward or downward trend or seasonal effects and, the mean or variance being 
consistent over time. Violation of stationarity in climate extremes can lead to wrong estimates of return levels derived from 
a given extreme value model [2]. There are many tests for stationarity in the statistical literature. In this study, we use the 
Kwiatkowski, Phillips, Schmidt and Shin test for stationarity introduced by in 1992. 

The KPSS test tests the null hypothesis that a time series is stationary (has no unit root) against the alternative that 
the series is non-stationary due to the presence of a unit root. Details of the test are summarised in [2]. Briefly, the test 
starts with a specification of the model 

x t — f?'D t + fx t + u t with fx t = y u t _ 1 + £ t and s t ~WN(0, ct 2 ), 

where fi is a vector of regression coefficients, D t contains deterministic components (constant or constant plus time trend), 
H t is a pure random walk with innovation variance er 2 , and fx t is a random error term of x t ,and may be heteroskedastic, 
£ t denotes the error term of /x t and is assumed to be a series of independent and identically distributed random variables 
with a mean of zero and constant variance er 2 . The quantity t is the time index for each function, with t — 1,2 ..., T. 

The test statistic for the KPSS test is a score statistic for testing er 2 = 0 against the alternative that er 2 > 0 and is 

— §2 I 

defined by K = r2 / where S 2 = 2y=i Uj, Uj is the residual of the regression of X t on D t and is a consistent 

estimate of the long-run variance of u t using u t . And, since the test is one-sided right-tailed, the null hypothesis of 
stationarity is rejected at the 100a% level if the KPSS test statistic K is greater than the 100(1 — a)% quantile from the 
appropriate asymptotic distribution. 

Probability Distributions 

Using the generalized Q-Q plot as a guide for determining the maximum domain of attraction of the underlying 
distribution for the annual maximum rainfall data for Botswana, for the period 1901-2012, we consider two probability 
distributions in the Gumbel domain (gamma and lognormal) and one in the Frechet domain (BurrXII). The justification for 
this choice is discussed under the results section of this paper.A description of each of the probability distributions under 
consideration is contained in Table 1. 


Table 1: Descriptions of Probability Distributions used in the Study 


Distribution 

Probability Density Function 

Range 

Parameters 

Gamma 

/(*) = ^«r(a) % “ lexp( * //?) 

0 < X < oo 

a > 0 is a continuous shape 
parameter; /? > 0 is a 
continuous scale parameter 

Lognormal 

fix) = /—exp ( (log(x) n) 2 /2a 2 ) 

xa^zn 

0 < X < oo 

— oo < ix < oo is the mean 
and a > 0 are the mean and 
standard deviation of the 
natural logarithm of X 

BurrXII 

(XT f X \ a ~ 1 ( /'X\“') -(T+1 ' ) 

K-)} 

0 < X < oo 

(3 > 0 is the scale parameter, 
oc> 0 is a continuous shape 
parameter and r > 0 is a 
continuous shape parameter 
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The most popular maximum likelihood estimation method is used to calculate the parameter estimates for the 
models under study using various packages of the R language. 

Goodness-of-Fit and Model Selection 


To assess the fit of the models under consideration, two diagnostic plots, namely, the empirical cumulative 
distribution function (CDF) and the Quantile plots, and two more formal criteria, the Anderson-Darling (AD) test and 
Bayesian information criterion (BIC) are used. 


The AD test tests the null hypothesis that a sample of size n comes from a population with a specified distribution. 
It is based on the discrepancy between the empirical cumulative distribution function (ECDF), F n (x), and the theoretical 
cumulative distribution function (CDF), F(x ), of a specified distribution. This test gives more weight to the tails of the 
distribution than the more popular Kolmogorov-Smirnov test. As explained -by [19], in the case of fitting maximum 
values, if x^,x^, ■■■,x^ n - ) denote an increasingly ordered sample of independent observations from a population with 
distribution function F(x), the Anderson-Darling statistic is calculated using Equation (2.4). 


A 2 


r)^°,g(r z n+i-D 
' n 


■Z?=i z i / n 


(2.4] 


where Zj = F 0 {x^), the fitted distribution function evaluated at the ith element of the ordered sample of size n. The best-fit 
model (the model that gives the best fit to the data) is the one with the smallest value of the AD test statistic under 
consideration. 


The BIC is an information theoretic criterion for model selection among a finite set of models with different 
numbers of parameters. It is given by the formula in Equation (2.5). 

BIC = -2lnf(x\§) + dlnn (2.5] 

where f{x |$) is the score of the maximum likelihood function for the estimated model, n is the sample size and d 
is the total number of parameters. The model with the smallest value of the BIC is the most plausible. 

RESULTS AND DISCUSSIONS 

Figure 1 shows a time series plot and histogram of the maximum annual average rainfall for Botswana 
over the years 1901-2012. The histogram shows that the data have moderate positive skewness. 




O SO -too 150 200 250 

Annual maximum rail nrsui<mm) 


Figure 1: Time Series Plot (Left Panel) and Histogram (Right Panel) of 
Annual Maximum Rainfall for Botswana, 1901-2012 
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Descriptive Statistics 


Table 2: Descriptive Statistics 


N 

Mean 

Stdev 

Ql 

Median 

Q3 

Min 

Max 

Skew 

Kurtosis 

Se 

112 

112.12 

12 36 

83.45 

36.11 

134.18 

51.02 

233.58 

0.76 

0.55 

3.41 


The results contained in Table 1, indicate that during the years 1901 to 2012, Botswana received the highest 
annual maximum rainfall of 233.60 mm and the lowest of 51.02 mm. The value of excess kurtosis of y 2 = 0.55 is greater 
than that of the normal distribution (which is 0), shows that the data are heavy-tailed relative to the normal distribution. 
And, the coefficient of skewness of y 1 = 0.76 indicates that the data are positively skewed, as shown by the histogram in 
Figure 1. Further to the summary statistics contained in Table 1, a Jarque-Bera test for normality produces a chi-squared 
value of x 2 — 12.939 on 2 degrees of freedom, with a p — value — 0.00155, indicating that the skewness and kurtosis 
for these data significantly differ from their expected values under normality. 

Testing for Stationarity 

The KPSS test gives aKPSSLevel — 0.1285, Truncationlagparameter = 4, with ap — value = 0.1 
(Warning message: In kpss.test (y): p-value greater than printed p-value). And, since this p-value is significantly large, we 
fail to reject the null hypothesis of stationarity and conclude that the Botswana annual maximum rainfall data for the stated 
period are stationary. 

Selecting Maximum Domain of Attraction 

The generalized Q-Q plot is shown in Figure 2. The plot shows a fairly constant behaviour or trend in the large 
observations, suggesting the underlying distribution of the data is most probably in the Gumbel domain. This pattern is 
then followed by an increasing behavior in the largest 4 quantiles, which suggests that a distribution in the Frechet domain 
might give a reasonably good fit to the annual maximum rainfall for Botswana data. As a result, we focus on probability 
distributions in the Gumbel domain and, consider two distributions that are widely used for modelling rainfall data: the 
gamma and lognormal. For comparison, and still based on the behavior of the generalized Q-Q plot in the largest 
observations, we fit the three-parameter Burr XII, which belongs in the Frechet domain, to the data. 



Quantiles of standard exponential 

Figure 2: Generalized Q-Q Plot for the Annual Maximum Rainfall for Botswana, 1901-2012 


www.iaset.us 


editor@iaset. us 
















8 


Wilson Moseki Thupeng 


Parameters Estimates 

Table 3 contains maximum likelihood estimates (MLEs) of the parameters of the three fitted models for the 
annual maximum rainfall for Botswana for the stated period. The standard error of each estimate is given in brackets. The 
standard errors for the gamma and lognormal parameter estimates are generally small, indicating small estimation variance. 
This is in sharp contrast to the standard errors of the estimates of the BurrXII, which are relatively large, especially in the 
case of the scale parameter/?. This suggests that the BurrXII is the least preferred model of the three models considered for 
the Botswana annual maximum rainfall data recorded during the period 1901-2012. The next section discusses the 
goodness-of-fit of the three models. 


Table 3: MLE Parameter Estimates of the Fitted Models 


Distribution 

Estimate (Standard Error) 

Gamma(2P) 

a = 10.13566663 (1.33061037) 

(3 = 0.09040273 (0.01216516) 

Lognormal (2P) 

p = meanlnx = 4.6694692 (0.03001591) 
ct = Sdlnx = 0.3176586 (0.02122351) 

BurrXII (3P) 

a = shapel = 1.459122 (0.7112521) 
x = shape2 = 4.835385 (0.7298173) 

P = Scale = 119.664023 (18.5060675) 


Goodness-of-Fit and Model Selection 

The Empirical versus theoretical distribution functions (CDF) plot (left) and the Q-Q plot (right) for the gamma, 
lognormal and Burr XII models are shown in Figure 3 


6 




data 


Theoretical quantiles 


Figure 1: Empirical Versus Theoretical Distribution Functions (Left) and Q-Q Plots for the Models 

From the CDF plot, it appears that the three models all fit the data very well. However, the Q-Q plot makes things 
a bit clearer and indicates that the gamma and the lognormal are very competitive best-fit models, with the BurrXII ranking 
a distant third. 


The findings based on the diagnostic plots in figure 3 are confirmed by the results of the AD test and the BIC 
contained in Table 2. Both the AD test and the BIC agree on the ranking of the three models: the lognormal distribution 
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gives the best fit to the data followed closely by the gamma whilst the BurrXII distribution comes a distant third. It must be 
noted, however, that based on both the diagnostic plots in Figure 3 and goodness-of-fit criteria in Table 4, the BurrXII also 
gives a reasonably good fit to the annual maximum rainfall data under consideration. Overall, these results are to be 
expected given the tail behaviour of these data as indicated by the generalized Q-Q plot in Section 3.3 of this paper. 


Table 4: The A-D Test Statistics and BIC Scores for the Fitted Models 


Goodness-of-Fit Statistic/Criterion 

BurrXII 

Gamma 

Lognormal 

Anderson-Darling 

0.54129838 

0.42747640 

0.40346494 

Rank 

3 

2 

1 

Bayesian Information Criterion 

1124.860 

1117.472 

1116.362 

Rank 

3 

2 

1 


CONCLUSIONS 

In this paper, annual maximum rainfall data for Botswana for the period January 1901 to December 2012 are 
modelled using the gamma, lognormal and BurrXII probability distribution. Application of the Kwiatkowski-Philips- 
Schmidt-Shin test of stationarity on the data shows that the time series is stationary. The maximum likelihood method is 
used to estimate the parameters of each distribution under study. The CDF probability plot and Q-Q plot show that the 
three models are highly competitive in terms of goodness-of-fit. Moreover, both the AD test and the BIC agree on the 
ranking of the three models: the lognormal distribution gives the best fit to the data followed closely by the gamma whilst 
the BurrXII distribution comes a distant third. 
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